﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Runtime.CompilerServices;
using System.Text;
using RayDen.Library.Core.Primitives;

namespace RayDen.Library.Core
{
    using Library.Core.Primitives;
    using RNG = FastRandom;

    /// <summary>
    /// Monte-Carlo routines
    /// </summary>
    public class MC
    {
        const float OneMinusEpsilon = 0.9999999403953552f;


        // Sampling Function Definitions
        public static void StratifiedSample1D(float[] samp, int nSamples, RNG rng, bool jitter)
        {
            float invTot = 1f / nSamples;
            for (int i = 0; i < nSamples; ++i)
            {
                float delta = jitter ? rng.NextFloat() : 0.5f;
                samp[i] = Math.Min((i + delta) * invTot, OneMinusEpsilon);
            }
        }

        public static void StratifiedSample1D(double[] samp, int nSamples, RNG rng, bool jitter)
        {
            double invTot = 1.0 / nSamples;
            for (int i = 0; i < nSamples; ++i)
            {
                var delta = jitter ? rng.NextDouble() : 0.5;
                samp[i] = Math.Min((i + delta) * invTot, OneMinusEpsilon);
            }
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static void StratifiedSample2D(float[] samp, int nx, int ny, RNG rng, bool jitter)
        {
            float dx = 1f / nx, dy = 1f / ny;
            int samps = 0;
            for (int y = 0; y < ny; ++y)
                for (int x = 0; x < nx; ++x)
                {
                    float jx = jitter ? rng.NextFloat() : 0.5f;
                    float jy = jitter ? rng.NextFloat() : 0.5f;
                    samp[samps++] = Math.Min((x + jx) * dx, OneMinusEpsilon);
                    samp[samps++] = Math.Min((y + jy) * dy, OneMinusEpsilon);
                }
        }

        static FastRandom Randomizer = new FastRandom();
        public static void RejectionSampleDisk(ref float x, ref float y)
        {
            float sx, sy;
            do
            {
                sx = 1f - 2f * (float)Randomizer.NextDouble();
                sy = 1f - 2f * (float)Randomizer.NextDouble();
            } while (sx * sx + sy * sy > 1f);
            x = sx;
            y = sy;
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static Vector UniformSampleHemisphere(float u1, float u2)
        {
            float z = u1;
            float r = MathLab.Sqrt(Math.Max(0f, 1f - z * z));
            float phi = 2f * MathLab.M_PI * u2;
            float x = r * (float)Math.Cos(phi);
            float y = r * (float)Math.Sin(phi);
            return new Vector(x, y, z);
        }
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static float UniformHemispherePdf(float theta, float phi)
        {
            return MathLab.INVTWOPI;
        }
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static Vector UniformSampleSphere(float u1, float u2)
        {
            float z = 1f - 2f * u1;
            float r = MathLab.Sqrt(Math.Max(0f, 1f - z * z));
            float phi = 2f * MathLab.M_PI * u2;
            float x = r * (float)Math.Cos(phi);
            float y = r * (float)Math.Sin(phi);

            return new Vector(x, y, z).Normalize();
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static float UniformSpherePdf()
        {
            return 1f / (4f * MathLab.M_PI);
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static void UniformSampleDisk(float u1, float u2, ref float x, ref float y)
        {
            float r = MathLab.Sqrt(u1);
            float theta = 2.0f * MathLab.M_PI * u2;
            x = r * (float)Math.Cos(theta);
            y = r * (float)Math.Sin(theta);
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static void ConcentricSampleDisk(float u1, float u2, ref float dx, ref float dy)
        {
            float r, theta;
            // Map uniform random numbers to $[-1,1]^2$
            float sx = 2 * u1 - 1;
            float sy = 2 * u2 - 1;
            // Map square to $(r,\theta)$
            // Handle degeneracy at the origin
            if (sx == 0.0 && sy == 0.0)
            {
                dx = 0.0f;
                dy = 0.0f;
                return;
            }
            if (sx >= -sy)
            {
                if (sx > sy)
                {
                    // Handle first region of disk
                    r = sx;
                    if (sy > 0.0)
                        theta = sy / r;
                    else
                        theta = 8.0f + sy / r;
                }
                else
                {
                    // Handle second region of disk
                    r = sy;
                    theta = 2.0f - sx / r;
                }
            }
            else
            {
                if (sx <= sy)
                {
                    // Handle third region of disk
                    r = -sx;
                    theta = 4.0f - sy / r;
                }
                else
                {
                    // Handle fourth region of disk
                    r = -sy;
                    theta = 6.0f + sx / r;
                }
            }
            theta *= MathLab.M_PI / 4f;
            dx = r * (float)Math.Cos(theta);
            dy = r * (float)Math.Sin(theta);
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static void UniformSampleTriangle(float u1, float u2, ref float u, ref float v)
        {
            float su1 = MathLab.Sqrt(u1);
            u = 1f - su1;
            v = u2 * su1;
        }
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static float UniformConePdf(float cosTheta)
        {
            return 1f / (2f * MathLab.M_PI * (1f - cosTheta));
        }



        //public static Vector UniformSampleCone(float u1, float u2, float costhetaMath) {
        //    float costheta = Lerp(u1, costhetaMath, 1f);
        //    float sintheta = MathLab.Sqrt(1f - costheta*costheta);
        //    float phi = u2 * 2f * MathLab.PI;
        //    return Vector((float)Math.Cos(phi) * sintheta,
        //                  (float)Math.Sin(phi) * sintheta,
        //                  costheta);
        //}

        public static void LatLongMappingMap(float s, float t, out Vector wh, out float pdf)
        {
            float theta = t * MathLab.M_PI;
            float sinTheta = (float)Math.Sin(theta);

            float phi = s * 2f * MathLab.M_PI;
            wh = MathLab.SphericalDirection(sinTheta, (float)Math.Cos(theta), phi);
            pdf = MathLab.INVTWOPI * MathLab.INVPI / sinTheta;
        }

        public static Vector CosineSampleHemisphere(float u0, float u1)
        {
            var ret = new Vector();
            ConcentricSampleDisk(u0, u1, ref ret.x, ref  ret.y);
            ret.z = MathLab.Sqrt(Math.Max(0f, 1f - ret.x * ret.x - ret.y * ret.y));
            return ret;
        }


        public static Vector CosineSampleHemisphere(float u0, float u1, out float pdfW)
        {
            var ret = new Vector();
            ConcentricSampleDisk(u0, u1, ref ret.x, ref  ret.y);
            ret.z = MathLab.Sqrt(Math.Max(0f, 1f - ret.x * ret.x - ret.y * ret.y));
            pdfW = ret.z * MathLab.INVPI; 
            return ret;
        }

        public static float Mutate(RNG RandomProvider, float value, float t1 = 1f/1024f, float t2 = 1f/64f)
        {
            var returnValue = value;
            var s1 = t1;
            var s2 = t2;
            var dv = s2 * Math.Exp(-Math.Log(s2 / s1) * RandomProvider.NextDouble());
            if (RandomProvider.NextDouble() < 0.5)
            {
                returnValue += (float)dv;
                if (returnValue > 1f)
                    returnValue -= 1f;
            }
            else
            {
                returnValue -= (float)dv;
                if (returnValue < 0f)
                    returnValue += 1f;
            }
            return returnValue;
        }

        public static float Mutate(float u0, float u1, float value, float t1 = 1f/1024f, float t2 = 1f/64f)
        {
            var returnValue = value;
            var s1 = t1;
            var s2 = t2;
            var dv = s2 * Math.Exp(-Math.Log(s2 / s1) * u0);
            if (u1 < 0.5)
            {
                returnValue += (float)dv;
                if (returnValue > 1f)
                    returnValue -= 1f;
            }
            else
            {
                returnValue -= (float)dv;
                if (returnValue < 0f)
                    returnValue += 1f;
            }
            return returnValue;
        }
        public static double Mutate(double u0, float u1, float value, float t1 = 1f/1024f, float t2 = 1f/64f)
        {
            var returnValue = value;
            var s1 = t1;
            var s2 = t2;
            var dv = s2 * Math.Exp(-Math.Log(s2 / s1) * u0);
            if (u1 < 0.5)
            {
                returnValue += (float)dv;
                if (returnValue > 1f)
                    returnValue -= 1f;
            }
            else
            {
                returnValue -= (float)dv;
                if (returnValue < 0f)
                    returnValue += 1f;
            }
            return returnValue;
        }
        public static void ComputeStep1dCDF(float[] f, int nSteps, float[] cdf, out float c)
        {
            cdf[0] = 0f;
            for (int i = 1; i < nSteps + 1; ++i)
                cdf[i] = cdf[i - 1] + f[i - 1] / nSteps;
            c = cdf[nSteps];
            for (int i = 1; i < nSteps + 1; ++i)
                cdf[i] /= c;
        }


        public static float BalanceHeuristic(int nf, float fPdf, int ng, float gPdf)
        {
            return (nf * fPdf) / (nf * fPdf + ng * gPdf);
        }
        public static float PowerHeuristic(int nf, float fPdf, int ng, float gPdf)
        {
            float f = nf * fPdf, g = ng * gPdf;
            return (f * f) / (f * f + g * g);
        }

        public static void Shuffle<T>(T[] samp, int count, int dims, FastRandom rng)
        {
            for (var i = 0; i < count; ++i)
            {
                var other = i + (rng.NextInt() % (count - i));
                for (var j = 0; j < dims; ++j)
                    Swap(ref samp[dims * i + j], ref samp[dims * other + j]);
            }
        }

        public static double RadicalInverse(int n, int _base)
        {
            double val = 0;
            double invBase = 1.0 / _base, invBi = invBase;
            while (n > 0)
            {
                // Compute next digit of radical inverse
                int d_i = (n % _base);
                val += d_i * invBi;
                n = (int)(n * invBase);
                invBi *= invBase;
            }
            return val;
        }

        
        public static void Sample02(uint n, uint[] scramble,float[] sample) {
            sample[0] = VanDerCorput(n, scramble[0]);
            sample[1] = Sobol2(n, scramble[1]);
        }
        public static void Sample02(uint n, uint scramble0, uint scramble1, out float sample0, out float sample1)
        {
            sample0 = VanDerCorput(n, scramble0);
            sample1 = Sobol2(n, scramble1);
        }

        public static  float VanDerCorput(uint n, uint scramble) {
            // Reverse bits of _n_
            n = (n << 16) | (n >> 16);
            n = ((n & 0x00ff00ff) << 8) | ((n & 0xff00ff00) >> 8);
            n = ((n & 0x0f0f0f0f) << 4) | ((n & 0xf0f0f0f0) >> 4);
            n = ((n & 0x33333333) << 2) | ((n & 0xcccccccc) >> 2);
            n = ((n & 0x55555555) << 1) | ((n & 0xaaaaaaaa) >> 1);
            n ^= scramble;
            return Math.Min(((n>>8) & 0xffffff) / (float)(1 << 24), OneMinusEpsilon);
        }


        public static float Sobol2(uint n, uint scramble) {
            for (uint v = 1u << 31; n != 0; n >>= 1, v ^= v >> 1)
                if ((n & 0x1u)== 1) scramble ^= v;
            return Math.Min(((scramble>>8) & 0xffffff) / (float)(1 << 24), OneMinusEpsilon);
        }


        public static float LarcherPillichshammer2(uint n, uint scramble) {
            for (uint v = 1u << 31; n != 0; n >>= 1, v |= v >> 1)
                if ((n & 0x1)==1) scramble ^= v;
            return Math.Min(((scramble>>8) & 0xffffff) / (float)(1 << 24), OneMinusEpsilon);
        }

        private static void Swap<TItem>(ref TItem t1, ref TItem t2)
        {
            TItem swp = t1;
            t1 = t2;
            t2 = swp;
        }

        public static float Mutate(double p1, float p2, float p3)
        {
            return Mutate(p2, p3,(float) p1);

        }
    }
}
